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Abstract: We present a next-to- leading order calculation of Higgs boson production plus 
one and two jets via gluon fusion interfaced to shower Monte Carlo programs, implemented 
according to the POWHEG method. For this implementation we have used a new interface 
of the POWHEG BOX with MadGraph4, that generates the codes for generic Born and real 
processes automatically. The virtual corrections have been taken from the MCFM code. We 
carry out a simple phenomenological study of our generators, comparing them among each 
other and with fixed next-to-leading order results. 
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1. Introduction 

The search for the Standard Model Higgs boson is entering the endgame phase. Since the 
start of the LHC running, the allowed mass range for the Standard Model Higgs boson has 
been greatly curtailed. Furthermore there are tantalizing, but inconclusive, hints in the 
remaining low mass region [1, 2]. Should these hints be confirmed by further data, it will 
be a matter of some urgency to examine the properties of the new state (or states). In the 
low mass region, the Standard Model Higgs boson is predicted to have about ten decay 
modes with branching fractions greater than one per mille, so there will be a number of 
channels to be studied. 

In order to extract the Higgs boson couplings to fermions and gauge bosons from these 
different channels, detailed information about the production rates will be required [3, 4, 5]. 
In this context, the Higgs boson + 2 jet process enters principally as an irreducible back- 
ground for the vector boson fusion (VBF) process [6, 7, 8]. This process is of particular 
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importance in the extraction of the HWW and HZZ couphngs and in confirming if tlie 
detected scalar particle is the Higgs boson responsible for electroweak symmetry breaking. 
In order to distinguish VBF Higgs boson signal from backgrounds, stringent cuts are re- 
quired on the Higgs boson decay products as well as on the two forward quark jets which 
are characteristic for VBF. Information on the CP properties of the Higgs boson can be 
extracted by studying the azimuthal distributions of the two hardest jets [4, 9], in events 
with a large (> 3 units) rapidity separation between them. In fact, the signature of the 
VBF Higgs boson production are two jets well separated in rapidity. In addition, we expect 
very low hadronic activity between the two hardest jets, due to the exchange of colour- 
less vector bosons in the t channel, in contrast to what is expected for the gluon-fusion 
production mechanism. For the gluon-fusion processes, the higher-order corrections are 
significant [10, 11, 12, 13], and a detailed understanding of the structure of the radiation 
pattern is needed in order to access the efficiency of the central-jet veto [14, 15]. 

The Higgs boson search channels are often subdivided according to the number of 
associated jets, for a number of reasons. First, a differing number of associated jets can 
imply a different source of background. For example, in the search for a Higgs boson 
decaying to W^W~ , the backgrounds have different origins in the differing jet bins, so it 
proves advantageous to analyze the different jet multiplicities separately. Second, in the 
two-jet bin, the new VBF production channel will come into play. For these reasons it is 
pressing to provide NLO predictions for gluon-fusion initiated Higgs boson -|- jets processes 
in the POWHEG formalism, so that experimentalists can incorporate the best theoretical 
information into their analyses. 

Phenomenological studies for the production of a Higgs boson in association with two 
jets, at the parton level and at fixed order, have been available for a while. In refs. [14, 15], 
the calculation was performed at leading order, with the gluons coupling to the Higgs 
boson via a top-quark loop, retaining the exact top-quark mass dependence (mj) in the 
whole calculation. In that paper it was shown that the large-mj limit provides an excellent 
approximation to the full rrit dependence when the Higgs boson mass, uih, is small com- 
pared to the top-pair threshold. The large-mj limit was found to break down for > mt 
and when jet transverse momenta become large {p{, > nit). However, large dijet invariant 
masses do not invalidate the nit ^ oo limit, as long as the Higgs boson mass and the 
jet transverse momenta are small enough, less than the top-quark mass in practice. This 
observation opened the possibility to compute the NLO corrections to Higgs boson -|- 2 jet 
production in gluon fusion, in the large-mt limit, starting with an effective coupling of the 
Higgs boson to the gluon field [12, 13] 

In this work we have implemented Higgs boson plus one- and two-jet production in 
gluon fusion, in the large top-quark mass limit, using the POWHEG method. While for the first 
process a matched NLO-|-shower calculation has already appeared in the literature [16], the 
second one has never been performed before. We have built our generators using the POWHEG 
BOX framework [17], with the virtual corrections taken from the MCFM program [18], and the 
Born , colour-correlated Born, spin-correlated Born and real contributions computed using 
a new POWHEG BOX interface to the MadGraph4 [19, 20] program, which is fully generic and 
can be applied to any process. Thus, the aim of this paper is twofold: 
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• To present a new generic interface of MadGraph4 to the POWHEG BOX that allows for the 
automatic generation of the code for the Born amplitude, the Born colour- and spin- 
correlated contributions, for the real-radiation amplitude and for the Born colour 
flow in the leading-colour approximation. These ingredients are all needed in the 
implementation of a new process in the POWHEG BOX. This interface can be used for 
any process that can be generated with MadGraph4. Thanks to this interface, in order 
to construct a POWHEG BOX generator, one needs only to provide the Born phase space 
and the virtual corrections. 

• To illustrate the Higgs plus one and two jet generators, by comparing their outputs 
to the corresponding NLO results and amongst themselves. 

The paper is organized as follows: in sec. 2, we present the new interface of the POWHEG BOX 
to MadGraph4. In sec. 3, we discuss the approximations used in our calculation and how 
wc extract the virtual corrections from the MCFM code. In sec. 4 we give more details about 
the phase-space generator for Higgs boson plus one and two jets, and how we deal with 
the divergences present at the Born level. After briefly recalling a few crucial POWHEG BOX 
parameter settings, we discuss some phenomenology in sec. 5. In particular we compare 
the NLO differential cross sections with the results obtained with the POWHEG BOX, at 
several levels of approximation. Finally, in sec. 6, we summarize our findings. A few 
technical details of the MadGraph4 interface and of the PYTHIA setup are collected in the 
two appendices. 

2. The interface to MadGraph4 

The MadGraph4 package [19, 20] generates squared tree-level matrix elements in an auto- 
matic way. It is based on Feynman diagrams, using HELAS routines [21] as building blocks 
which are subsequently combined into colour-ordered amplitudes and written in a FOR- 
TRAN code. Furthermore, the MadDipole [22, 23, 24] extension of MadGraph4 and the 
MadFKS tool [25] can generate the spin- and colour-correlated Born squared amplitudes. 
The MadGraph4 package can thus construct all the ingredients needed to interface a NLO 
calculation to the POWHEG BOX (see the complete list in the introduction of ref. [17]), with 
the exception of the Born phase space and the virtual matrix elements. This makes the 
MadGraph4 framework an ideal environment to complement the POWHEG BOX in the creation 
of all tree-level squared matrix elements and in this section we describe a newly developed 
interface between these two codes. 

For a given process, the routines and parameters that are generated by MadGraph4 for 
the POWHEG BOX are (see ref. [17] for more details on the notation and conventions used): 

• the multiplicity of the Born and real-emission processes, nlegborn and nlegreal, 
respectively, defined in the nlegborn. h include file; 

• the list of Born and real-emission flavour structures 
flst_born(k=l , nlegborn, j=l : f lst_nborn) , 
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f lst_real (k=l ,nlegreal , j=l : f lst_nreal) , 

respectively, to be initialised by a call to the init_processes routine; 

• the routine 

setborn(p(0 : 3, 1 :nlegborn) ,bf lav(l :nlegborn) ,born, 

bornjkd :nlegborn, 1 :nlegborn) ,bmunu(0 : 3 , : 3, 1 :nlegborn) ) 
which computes, for a given set of four-momenta p and flavour structure bf lav, 
the Born squared matrix element born, the colour-correlated, bornjk, and the spin- 
correlated one, bmunu; 

• the routine that computes the real emission squared matrix elements amp2 (stripped 
of a factor as/(27r)), for a given momentum configuration p and flavour structure 
rf lav 

setreal(p(0:3,nlegreal) ,rf lav(l :nlegreal) ,aiiip2); 

• a colour flow assignment to the Born squared matrix elements, in the leading colour 
approximation, needed by the showering programs. The colour is assigned on statisti- 
cal grounds, based on information that is cached during the call to the squared matrix 
elements, by the routine borncolour_lh. In addition to this routine, the interface 
also provides a similar one for the real squared matrix elements, realcolour_lh. In 
the present version of the POWHEG BOX this routine is not used; instead a different 
method is employed to assign the colour to the real-radiation matrix elements, as 
discussed in sec. 8 of ref. [17]; 

• an interface to the Lcs Houches parameter input card to specify the physics-model 
parameters through the routine init_couplings. 

In contrast to the default way of generating amplitudes with MadGraph4, all the squared 
matrix elements for the various flavour structures are written in the same directory, making 
sure that the flies and routines have different names. This allows one to compile the code 
into a single library that contains all the matrix elements for all the flavour structures. 
An interface is provided that concatenates the flavour vectors bf lav or rf lav into strings, 
which are used as unique identifiers of each of the squared matrix elements. 

More technical details on the use of the MadGraph4 interface to the POWHEG BOX can 
be found in appendix A. 

3. Higgs boson production in gluon fusion 

The amplitudes for a Higgs boson in association with three, four or five partons, needed 
for the Higgs boson plus one and two jet cross sections at NLO, are calculated using an 
effective Lagrangian to express the coupling of the gluons to the Higgs field [26] 

C^ = ^HivG^,G^'\ (3.1) 
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where the trace is over the colour degrees of freedom. At the order required in this paper, 
the coefficient C is given in the MS scheme by [27, 28] 

where v is the vacuum expectation value of the Higgs field, v = 246 GeV. 

We have used the automatic interface described in sec. 2 to generate the code for the 
Born amplitude, the Born colour- and spin-correlated contributions and for the real cross 
section. In appendix A we give some more details on the adopted procedure. The one-loop 
amplitudes for the Higgs boson plus three and four parton processes are extracted from 
MCFM as described below. 

3.1 Virtual cross sections 

The complete set of one-loop amplitudes for all Higgs boson -|- 4 parton processes have been 
available [29, 30, 31, 32, 33, 34, 35] for some time. These formulae have been implemented 
into MCFM [13]. 

The interface between MCFM and the POWHEG BOX is fairly straightforward. Since in the 
POWHEG BOX implementation the Born, real and subtraction terms are available indepen- 
dently, MCFM needs only to return the pure virtual contributions. The interface to MCFM 
transfers the electroweak parameters, scale and scheme choices. Once these are established, 
calls to the virtual routines of MCFM can be made on an event-by-event basis. 

Since the routines that fill the values of the electroweak parameters (and the other 
information) are generic, this interface could be fairly easily extended to include other pro- 
cesses currently implemented in MCFM. However the normal MCFM routines are designed to 
return matrix elements that are summed over the flavour of the final-state partons. There- 
fore, in order to correctly interface to the POWHEG BOX, one must make small modifications 
to the MCFM code such that it returns the matrix elements for individual final-state flavour 
combinations. 

4. Implementation of the Hj and Hjj generators 

In the following, we use the notation H, Hj and Hjj to refer to Higgs boson generators 
that, at the Born level, describe the production of a Higgs boson plus zero, one and two 
partons. 

Using the automated MadGraph4 interface described in sec. 2, we have built the routines 

for the Born amplitude, the Born colour- and spin-correlated contributions and the real 
squared amplitudes, directly in the format required by the POWHEG BOX. The generation 
of the amplitudes is fast, taking only a few minutes for the Hjj process. The virtual 
corrections were extracted from the MCFM code, as described in sec. 3.1. At this point, the 
only missing ingredient for completing our generators is the Born phase space. 

4.1 Phase space for the Hj generator 

The phase space routine for the Hj Born process is trivial, since it is just a two-body 
phase space. It has been implemented with the possibility to activate an optional cut on 
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the transverse momentum of the Higgs boson, and with an effective importance sampUng 
of the small transverse-momentum region. The cut is necessary if one wants to generated 
unweighted events. Alternatively the cut can be set to zero, and the cross section is 
weighted with a suppression factor, that is a function of the underlying-Born kinematics. 
This factor is equal to 

where Pt is the Higgs boson transverse momentum in the underlying-Born configuration 
and Pq is set by the user ^ . Events are then generated uniformly in the cross section times 
F, and with a weight proportional to 1/F. The normalization of the weight is such that 
the cross section for events passing a set of cuts is given by the sum of all weights for the 
events passing the cut divided by the total number of generated events. 

4.2 Phase space for the Hjj generator 

The phase space for the Hjj generator was built using the same factorized phase space that 
POWHEG uses for the real kinematics. In other words, we treat the Hjj Born phase space 
as the real phase space for the Hj process with one extra emission. We have allowed the 
possibility of performing importance sampling in the regions where the Born cross section 
becomes singular. In more detail, labelling 1 and 2 as the two final-state light partons in 
the Hjj process, we write the Hjj phase space using the following identity 

N El N Eo N N 

d^HU = d^HU -1- rp , r:. + d^H12 -1- „ , „ + d^H12 -j- + d^H12 -j- , (4.2) 

"12 J^l + ^2 "12 ^1 + -C/2 «1 "2 

where Ei is the energy of the i-th parton in the center-of-mass frame, and 

1 1 1 ^ 



where di, d2 and di2 are phase-space functions that vanish respectively when parton 1 is 
collinear to the initial-state partons, parton 2 is coUinear to the initial-state partons, and 
partons 1 and 2 arc collinear to each other. Their precise form is given in eqs. (4.23)-(4.26) 
of ref. [17]. We then factorize each phase-space factor in (4.2) according to the formula 

(21) N El (12) N E2 

d^H12 = d^H2 d^lJ — + d^Hl d^'.J — 

"12 -C/1 + -£^2 "12 -C/1 + -C'2 

+ d^H2 d^fi I + d^Hi £ . (4.4) 

The subscripts i^l or H2 characterize the underlying Born, while the superscript in $rad 

specifies the radiation process. Thus, for example, d^H2 is the underlying-Born phase 

(21) 

space of the Higgs boson together with parton 2, and d^^ad ^^'^ radiation phase space 
corresponding to parton 2 emitting parton 1. The notation d^rad"* ^^rad"* ^^ans that 



^The value of po, set to 20 GeV for the simulation done in this paper, or the form of the suppression factor 
F in eq. (4.1), can be changed by the user by modifying the bom_suppression routine in the Bom_phsp.f 
file. 
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partons 1 or 2 are emitted by the initial-state partons. The factorization of the phase space 
into underlying Born and radiation phase space for both initial- and final-state radiation 
is the default one used in POWHEG, and is described in detail in sees. 5.1 and 5.2 of ref. [36]. 
The decomposition in eq. (4.4) is such that appropriate importance sampling is performed 
in all singular regions. In fact, the factors in each term damp all but one singular region, 
and the corresponding factorized phase space performs importance sampling precisely in 
that region. Ideally, the phase-space integration should be performed with an integrator 
that can sum over a discrete variable. The POWHEG BOX integrator [37] does not have this 
feature at the moment. Thus, we divided the range of one extra integration variable into 
four segments, mapping each segment to one of the phase space components. 

The phase space of eq. (4.4) is not the default one used in the Hjj generator. It is 
activated by setting the variable fullphsp in the powheg_input file. The default phase 
space is simply 

d^Hi2 = d^m , (4.5) 

with no importance sampling at all. We have in fact observed that the loss of efficiency 
due to the increased number of calls to the matrix-element routines overwhelms any benefit 
arising from the improved importance sampling. 

As for the case of the Hj generator, the Hjj generator also includes the implementation 
of a Born suppression factor for the suppression of the singularities of the underlying-Born 
amplitude, in the born_suppression routine. It has the form 

F=( ^-M V (46) 

where po is a parameter that characterizes the minimum jet energy where some accuracy 
is required, Pti and the transverse momenta (with respect to the beam axis) of the 

two final-state partons, and 

p2^, = 2(l-cosM-^^, (4.7) 

that can be interpreted as the transverse momentum of parton 1 with respect to parton 2 
or vice-versa, depending upon which is the softest. This suppression factor can also play 
the role of a generation cut, if po is chosen small enough. 

4.3 The POWHEG BOX parameter setting 

We have turned on the bornzerodamp flag by default in both the Hj and the Hjj gen- 
erators. The purpose of this flag is explained in ref. [17]. This results in a considerable 
speedup of the Hjj code. However, no appreciable differences were observed in the results 
obtained without the bornzerodamp option. 

The separation of the singular regions is controlled in the POWHEG BOX by the param- 
eters par_diexp and par_dijexp that are set to 1 by default (see ref. [17], section 4.7). 
For the results presented here we have chosen to set them to 2, which leads to slightly 
better stability for the set of distributions that we have considered. We note that these 
parameters were also set to 2 in the POWHEG BOX dijet generator [38]. Higher values of these 
parameters correspond to sharper separation of singular regions. 
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5. Phenomenology 

In this section, we present a phenomenological study of our generators. This study does 
not aim to produce results with reahstic experimental cuts, but rather to compare the 
predictions of the generators for some key observables. In particular, wc compare the 
generator results to the fixed NLO result and with the POWHEG output at the level of the 
generation of the hardest radiation (i.e. before interfacing the result to a parton shower 
program), after the shower with no hadronization effects, and after hadronization. 

We present results for the LHC running at 7 TeV, computed using the CTEQ6M 
parton distribution function (pdf) set [39]. The same calculations can easily be performed 
with any other available set [40, 41], but a study of pdf effects is beyond the scope of the 
present paper. Jets are reconstructed using the anti-Zc^ jet algorithm [42], with R = 0.5 
and default recombination scheme. Cuts of 20, 50 and 100 GeV on the final-state jets are 
considered. 

The factorization- and renormalization-scale choice deserves a more detailed discussion. 

In view of the large NLO corrections for these processes, the scale dependence is quite large. 
This is also a consequence of the fact that Higgs boson production in gluon fusion starts at 
order a^, so that the Hj and the Hjj processes are of order a| and a|, respectively. While 
for H production the natural scale choice is of the order of the Higgs boson mass, in the case 
of Hj production one may have a two-scale problem, if the transverse momentum of the jet 
is much smaller or much larger than the Higgs boson mass. In addition, besides the Higgs 
boson mass, one may also consider the Higgs boson transverse momentum, which may be 
more appropriate for very small or very large Higgs boson transverse momentum p", or the 
Higgs boson transverse mass = ^ rn^ + (p^)^, which may be more appropriate for large 
Higgs boson transverse momentum. For the Hjj case, there are even more possibilities. 
Although a full study of scale dependence would be very valuable, we will not perform 
it in the present work. In the Hj case we will limit ourselves to the fixed Higgs boson 
mass scale choice and to transverse momentum of the Higgs boson in the underlying-Born 
configuration. In the Hjj case, we will consider the fixed Higgs boson mass scale choice 
and to the Ht scale choice at the underlying-Born level, where 



and are the final-state parton transverse momenta in the underlying-Born kinematics. 
All results shown in the next sections have been computed with uih = 120 GeV and 
T„ = 0.00575 GeV. The Higgs boson momentum has been generated distributed according 
to a Breit-Wigner function, with fixed width. 

In the following, we will label "LHE" the results obtained at the level of the POWHEG 
first emission, "PY" the results obtained with the POWHEG-I-PYTHIA combination with the 
hadronization and underlying event switched off, and with "NLO" the fixed NLO results. 

For some observables, there is overlap among the various generators. For example, 
the Higgs boson transverse momentum, as well as the one-jet multiplicity and the trans- 
verse momentum of the leading jet, are described by both the H and the Hj generators. 
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The H generator gives a plausible description for these quantities also at small transverse 
momenta, since it provides a resummation of transverse-momentum logarithms that the 
Hj program does not provide. On the other hand, at large transverse momenta, the Hj 
generator has NLO accuracy, something that the H generator does not have. For the same 
quantities, the Hjj generator cannot be used, since it requires the presence of a second jet. 
The two-jet multiplicity, as well as the second hardest jet transverse-momentum distribu- 
tion, are provided by both the Hj and the Hjj programs, again with a different level of 
accuracy. A comparison of the generators in the regions where they overlap will be carried 
out in sec. 5.4. 



5.1 Results for Hj production 

We have generated a sample with 2M events at the fixed scale and 7.5M at the running 
scale. The event generation time is approximately 25 seconds for 1000 events on a typical 
CPU. 

We have simulated Hj production with two scale choices: fip = n^, = mn and the 
Pj. of the underlying-Born parton /ip = /^r = Pt^- All the following plots will come in 
pairs, with the left plot referring to /ip = /Ur = itih, and the right one to /ip = //r = 
p^^. We compare the fixed NLO results, the POWHEG hardest-emission results (LHE) and 
POWHEG+PYTHIA (PY) ones, where the hadronization and underlying-event effects in PYTHIA 
have been turned off (see appendix B). We show results for a few physical observables with 
three different cuts applied to the final-state jets: Pt cuts of 20, 50 and 100 GeV. 




Figure 1: One- and two-jet multiplicities for the POWHEG hardest-emission (LHE) results, the 
POWHEG events showered with PYTHIA (PY) and the fixed NLO result (NLO). The left plot is 
computed using /ip = /^r = , while in the right plot /ip — Mr = Pt^ ■ 



In fig. 1 we compare the jet multiplicity for one and two jets, at the NLO, LHE and 
PY levels, for jet cuts of 20, 50 and 100 GeV. We notice that the one-jet multiplicity is 
similar in the three results. The PYTHIA multiplicity is smaller than the LHE, which can 
be understood as a result of showering off the first jet. More marked differences can be 
seen in the cross sections with two jets: here the LHE result is clearly larger than the NLO 
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one, showing however a very different pattern as a function of the jet cut, depending upon 
the scale choice. We wih clarify the origin of these patterns later, when we discuss the 
transverse-momentum spectrum of the jets. We point out, however, that the NLO result 
for the two-jet multiplicity is of order a|, and thus has a marked scale dependence. In 
the left plot the scale is equal to the Higgs boson mass, whilst in the plot on the right the 
scale is of the order of the minimum jet p^- For a jet cut of 20 GeV these scales are widely 
separated, which explains the large difference in the NLO results. 




0.5 
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Figure 2: Transverse momentum of the Higgs boson. The notation is as in fig. 1. 



In fig. 2 we show the transverse momentum of the Higgs boson. The LHE, the PY 
and the NLO results are in remarkable agreement for this quantity, which is expected, 
due its inclusiveness. The Higgs boson rapidity distribution is plotted in fig. 3 while the 
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Figure 3: Rapidity distribution of the Higgs boson in the Hj process. The three sets of curves refer 
to the three cuts on the hardest jet, i.e. p'^ > 20, 50 and 100 GeV, from top to bottom respectively. 



hardest jet Pt in fig. 4. Both these quantities display a behaviour similar to the Higgs 
boson transverse momentum, again expected due to their inclusiveness. 

We now turn to more exclusive quantities, beginning with the transverse momentum 
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Figure 4: The hardest jet distribution in the Hj process. 



of the second jet, shown in fig. 5. The characteristic NLO behaviour, diverging at smah 
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Figure 5: Transverse-momentum distribution of the second hardest jet in the Hj process, with a 
20 GeV cut on the transverse momentum of the first jet. Inset: the low region. 



transverse momenta, is clearly visible in the inset. The LHE result displays, instead, the 
typical Sudakov damping at small p^i compensated by an increase at larger p^- The 
fully showered result is smaller, since further showering degrades the second jet transverse 
momentum. The two-jet multiplicity, displayed in fig. 1, is consistent with the value of the 
cross section around Pt = 20 GeV in fig. 2. From the upper inset in the figure, we see that, 
at this value of transverse momentum, the LHE cross section is above the NLO one for the 
fixed-scale choice (left plot), while it is near to it for the running-scale choice (right plot). 
For larger cuts, the LHE cross section remains above the NLO one in both cases, and 
the ratio grows in the running-scale case. The comparison of the fixed- and running-scale 
choice in the figure also has some subtle features that should be remarked upon. The 
running-scale choice yields smaller /larger scales for smaller /larger second jet Pt- The LO 
cross section is larger for smaller scale choices, so that the NLO corrections are smaller, and 
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vice- versa for larger scale choices. Thus, for small px, the NLO corrections are smaller than 
for large Pt- In POWHEG, one first generates the underlying-Born configuration, according 
to the cross section B, which is the total inclusive cross section at a fixed underlying- 
Born configuration. The radiation kinematics is then generated using a shower technique. 
As a result the whole distribution for the radiation is amplified by a X-factor equal to 
B/B [43, 44, 45]. In this case, due to the growth of the NLO corrections as a function of 
Pt, the amplification of the radiated-jet distribution due to the B/B K -factor increases as a 
function of the transverse momentum, a trend that is visible in the right figure. Conversely, 
no such effect is present for fixed scales. However, in this last case, one should recall that, 
in the LHE events, one power of as is effectively evaluated at the transverse momentum of 
the jet rather than at the fixed scale, thus yielding a decrease in the cross section which is 
also visible in the left plot. The larger value of the LHE cross section with respect to the 
NLO one is related to the large /C-factor, i.e. is due to the fact that the hardest radiation 
is amplified by a factor B/B. 
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Figure 6: Transverse-momentum distribution of the second hardest jet using the B ^ B option. 



Figure 6 is similar to fig. 5 except that the LHE result is obtained setting the bornonly 
flag to true in the POWHEG BOX input file. When this flag is true, the POWHEG BOX generator 
uses the Born cross section, rather than the NLO B function, to generate the underlying- 
Born configuration. We can see that, in this case, the LHE result is much closer to the NLO 
result, except for small transverse momenta, where the Sudakov damping becomes manifest. 
This results confirms that the enhancement of the transverse-momentum distribution of 
the second hardest jet is indeed due to the B/B iT-factor. 

In fig. 7 we show the second hardest jet transverse momentum for different cuts on the 
first jet Pt- From the figure it is clear that, when the second jet has transverse momentum 
above the first jet Pt cut, the leading jet is forced to have larger transverse momentum, 
and the cross section falls more rapidly. The comparison among the LHE, PY and NLO 
curves shows the typical pattern, with the NLO diverging at small transverse momenta, 
the LHE being instead suppressed in that region, but raising above the NLO result because 
of the B/B ivT-factor. 
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Figure 8: The relative transverse moraenta of the partons within the hardest jet p^ '^^ as defined 
in eq. (5.2). 



In fig. 8 we plot the relative transverse momenta of the partons within the hardest 
jet p^t'^^. This quantity is defined as follows: we perform a longitudinal boost to a frame 
where the rapidity of the first jet ji is zero. In this frame we compute 

pf-^ = Y^\Pl2iM ^ (5.2) 

where kj is the momentum of the first jet, and pi the momenta of the partons clustered 
within the first jet. This quantity, when computed at the LHE level, displays a marked 
difference from the corresponding NLO result. As discussed in ref. [38], this can be easily 
understood if we remember that the LHE result is suppressed by a Sudakov form factor, 
which requires that no harder radiation has been emitted from either the initial- or the 
final-state partons. The large bin at p^''^^ = in the LHE result is due to events in which 
there is only one parton in the jet, most likely initial-state radiation events. On the other 
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hand, the hardest radiation provided by PYTHIA should restore a shape closer to the NLO 
result, although amplified by the B/ B factor. Further amplification of the showered result 
is due to the fact that the shower uses a running coupling evaluated at the scale of the 
radiation, while the NLO result uses higher scales, and to the presence of multiple emissions 
in the shower. It is clear that this distribution, being determined mostly by the shower 
program, is quite sensitive to the shower model and tuning, and to the interface between 
the shower and POWHEG BOX. 



5.2 Results for Hjj production 

We have generated a sample with 2.5M events both at fixed and running scales. The event 
generation time is approximately 25 minutes for 1000 events on a typical CPU. 

We have run the Hjj program for two scale choices, /ip = //r = and fip = fiB. = Ht- 
All the following plots will come in pairs, the left one referring to the first scale choice, and 
the right one referring to the second scale choice. 
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Figure 9: Two- and three-jet multiplicities in Hjj production. In the left plot the scale is chosen 
equal to the mass of the Higgs. In the right plot the scale is taken equal to 



We begin by showing in fig. 9 the two- and three-jet multiplicities in the Hjj process. 
We find considerable agreement of the LHE, PY and NLO output for both the two- and 
three-jet multiplicities for the left plot, which corresponds to the choice of scale /ip = 

= rriu. The right plot corresponds to the H^ choice of scale. In this plot the two-jet 
multiplicity agrees for all cuts, while the three-jet multiplicity displays marked differences 
between the LHE results and the NLO ones. For the 100 GeV transverse-momentum 
cut, the fully showered result and the NLO also differ. The reason for these differences 
is the following: with the scale choice the iT-factor is near one, and thus the B/B 
amplification that usually enhances the spectrum of the radiated jet (that in this case is 
the third jet) is not effective. The Ht scale is considerably larger than m^, especially with 
a large Pt cut, which implies a reduction of the Born cross section and an increase of the 
JC-factor. 
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Figure 10: Transverse-momentum distribution of the Higgs boson in a Hjj sample with a 20 GeV 
Pt cut on the two hardest jets. The left and right plots use respectively the and Ht scales. 



In fig. 10 we show the transverse-momentum distribution of the Higgs boson in a Hjj 
sample with a 20 GeV Pt cut on the two hardest jets. Good agreement is found between the 
LHE, PY and NLO curves, except for small Higgs transverse momenta, where differences of 
the order of 30% are found. This is not surprising, in view of the required presence of two 
relatively soft jets. Radiation off these jets is Sudakov suppressed in the LHE result with 
respect to the NLO one. Since this radiation would deplete the jets, the LHE result suffers 
less depletion, and is thus larger. On the other hand, the shower degrades the energy of 
the jets lowering the cross section, an effect visible in the PY result. For large Higgs boson 
Pt, at least one of the jets is forced to be hard, and thus these effects lose importance. 
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In fig. 11 we show the Higgs boson rapidity distribution in Hjj for pj. cuts of 20, 50 and 
100 GeV on the two hardest jets. Again, since this is an inclusive distribution, it displays 
good agreement between the LHE, PY and NLO results. The ratio is only displayed for 



the 20 GeV cut. The three predictions become even more consistent for larger px cuts, as 
one would expect. 
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Figure 12: Transverse-momentum distribution of the hardest jet in the Hjj process. The left and 
right plots use respectively the ttih and Ht scales. 
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Figure 13: Transverse-momentum distribution of the second hardest jet in the Hjj process. The 
left and right plots use respectively the itih and Ht scales. 



In figs. 12 and 13 we display the transverse-momentum distribution for the first and 
second jet. The first distribution bears some similarity to the Higgs transverse momentum, 
for inclusiveness reasons. Observe that the small momentum disagreement between the 
LHE and the NLO result observed for the Higgs transverse momentum case is less evident 
in the hardest jet plot, and even less so in the second jet spectrum. This is explained by 
the fact that the same point in the abscissa corresponds to an increasing hardness of the 
event in the Higgs boson, first jet and second jet pr distributions. 

Turning to less inclusive quantities, we show in fig. 14 the transverse-momentum dis- 
tribution of the third hardest jet. We recognize here the typical behaviour of the LHE 
results, with the damping of the small- momentum growth found in the NLO result, and 
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Figure 14: Transverse-momentum distribution of the third hardest jet in the Hjj process. The 
left and right plots use respectively the ttih and scales. 



with the typical increase due to the B/B factor at higher momenta. As already anticipated 
in the discussion on the jet multiplicity, we notice that the increase is modest for the uih 
choice of scale, but is very large for the Ht scale. 
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Figure 15: Transverse-momentum distribution of the third hardest jet in the Hjj process for the 
transverse momentum cut on the two hardest jets equal to 20, 50 and 100 GeV. 



In fig. 15 we display the third-jet transverse momentum using three different cuts on 
the Pt of the first and second jet. The pattern is similar to what was already observed 
in the Hj case. In the left plot, where the scale is chosen equal to m^, we see a better 
concordance of the LHE and NLO results as the transverse momentum increases. 

In figs. 16 and 17 we show the p^''^ distribution for the hardest and second hardest 
jet respectively. Here too the pattern is similar to what already observed for the Hj case. 
The LHE distribution is strongly suppressed, due to the fact that small p^'"-' values also 
imply that no initial-state radiation has taken place above that scale, and the distribution 
is dominated by the shower effects, that are not obtained with the same scale choice as 
the NLO result. We note again that the B/B X-factor plays a role here, especially for the 
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distribution of the hardest jet in the Hjj process. The left and right plots use 
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Figure 17: p^^''^ distribution of the second hardest jet in the Hjj process. The left and right 
plots use respectively the tuh and Ht scales. 



right plots, that have a larger value because of the use of the Hx scale. 
5.3 Hadronization effects and matching ambiguities 

In this section we focus upon two topics: the effects of hadronization and the ambiguities 
related to matching the LHE result to the shower. In figs. 18 and 19 we display the 
transverse-momentum distribution of the radiated jet for different cuts on the first jet pr, 
for the Hj and Hjj generators respectively. 

The curves labelled "PYscalup" are obtained with a non-default determination of the 
scalup parameter to be set in the Les Houches interface, which limits the hardness of 
the radiation of the following shower. While normally in POWHEG scalup is set to the 
hardest momentum of the radiation that POWHEG generates, in the "PYscalup" results, we 
determine it by finding the smallest transverse momentum of the LHE, which can be either 
the transverse momentum of any parton relative to the beam axis, or of any parton relative 
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Figure 18: Transverse-momentum distribution of the second hardest jet in the Hj process for the 
Pt cut on the hardest jet equal to 20, 50 and 100 GeV. 
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to any other parton, computed in the center-of-mass frame. Because of the way the real- 
radiation contribution is separated into singular contributions, the two choices differ only 
for subleading configurations, and one expects only a minor effect due to this change. The 
plots in both figures confirm this expectation. Similarly, we see that hadronization and 
underlying-event effects (indicated with HAD in the figures) have a sizable impact on the 
distributions only for small transverse momenta, as expected. 

5.4 Comparison betw^een the H, Hj and Hjj generators 

In this section, we compare a few distributions that are described by more than one available 
POWHEG BOX generator. This comparison can be considered as a first step in the direction 
of merging POWHEG BOX samples with increasing number of jets. 

We first consider a comparison between the Higgs boson transverse-momentum spec- 
trum obtained using the H and the Hj generators. The H generator describes this dis- 
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Figure 20: Comparison of the Higgs boson transverse-momentum distribution computed with the 
H , Hj and HqT generators. The LHE Hj results are shown for the niH and the Higgs underlying- 
Born p"^ scale choice. 



tribution in the whole Pt range, with the Sudakov region obtained with next-to-leading 
logarithmic accuracy, and with leading-order accuracy for the high-momentum region. On 
the other hand, the Hj generator describes this distribution with NLO accuracy, but only 
for transverse momenta above the Sudakov region. 

In these comparisons, the H sample was obtained following the recommendation of 
ref. [46], i.e. the parameter hfact was set to mjy/1.2 in the powheg. input file. With 
this setup, the resulting Higgs boson distribution is in remarkable agreement with the 
output of the HqT program [47, 48, 49, 50]. In addition, we apply a /C-factor of 1.32 to 
the H result, in order to match the HqT total NNLO cross section (cr^^° = 10.85 pb, 
^NNLo _ ]^4_35 p]-,^_ The results are displayed in fig. 20, which demonstrates the agreement 
between the results for the H sample and the output of the HqT program, for the Higgs 
boson Pt distribution [46]. This comparison is performed using a scale equal to the Higgs 
boson mass, since HqT accepts only a fixed renormalization and factorization scale. We also 
show two predictions for the Hj generator, using both the fixed and running-scale choices. 
The Hj generator predictions with a fixed scale are also in excellent agreement with HqT 
for large transverse momenta. This is not surprising, since, in this kinematic region, it is 
using the same ©(af) result as HqT. At low px, the Hj results begin to feel the lack of 
soft-gluon resummation effects that are included in HqT. The Hj result computed with the 
dynamical scale shows a substantially similar pattern, undershooting the HqT one by about 
25%. 

Turning to the comparison of the Hj and Hjj generators for the transverse-momentum 
distribution of the second hardest jet, we notice that the Hj generator computes this distri- 
bution down to small values of the transverse momentum, since it includes the appropriate 



-20- 





10-1 


> 




'Ge 


10-2 




10-3 




10-4 






b 


10-'^ 




10-^ 




1.5 




1 




0.5 





LHE Hjj 
LHE Hj 



> 20 GeV 



50 100 150 200 250 300 350 400 
p? [GeV] 



10-1 

10-2 

10-3 

10-1 

10-5 

10-^' 
1.5 
1 

0.5 





LHE Hjj ^ 
LHE Hj ' 



p'^ > 20 GeV 



50 100 150 200 250 300 350 400 
P? [GeV] 
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and Hjj generators. The left plot uses the m^, scale in both the Hj and Hjj generators. In the 
right plot the Hj generator uses the Higgs underlying-Born transverse-momentum scale, while the 
Hjj generator uses Ht- 



Sudakov effects. At large transverse momenta, however, it only has leading-order accuracy. 
The Hjj program has instead NLO accuracy at large transverse momenta, but does not 
include fully resummed Sudakov effects at small p-r- We do not have any higher- accuracy 
calculation for this distribution, as in the previous case. The results are displayed in fig. 21. 
Notice that, for the mu scale choice, the matching is very good, something that we expect 
since the X-factor is close to one in this case. Matching with the running scales leads, 
instead, to non-negligible differences, although a stability plateau is present roughly above 
60 GeV, provided that the transverse momentum is not too large. 

6. Conclusions 

In the present work we have developed generators for the gluon-fusion production of a 
Higgs boson in association with one or two jets, in the large top-quark mass limit. We have 
examined the output of the two generators, comparing them with fixed NLO calculations, 
among each other, and with the fully inclusive Higgs boson POWHEG generator. The features 
of the distributions we have examined are all well understood and reflect our expectations 
for a typical POWHEG generator. In general, we find that quantities inclusive in the radiated 
jet (i.e. the second jet in H + 1 jet and the third jet in H + 2 jets) are in good agreement 
with the NLO result, while distributions in the radiated jet reflect the features of Sudakov 
suppression and NLO enhancement that are typically found in NLO generators matched 
with a shower. We see no indication of problems related to the increase in complexity 
when going from the inclusive Higgs boson production to the associated production with 
one and two jets, other than an increase in the amount of computer time required for the 
calculation. 

The development of the Higgs boson production code has been achieved using a new 
interface to the MadGraph4 code, that has also been presented in this work. The use of 
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this interface has considerably simpUfied the construction of the generator. The interface 
is fuhy generic, and so we expect that the implementation of new processes will be greatly 
simplified with its use. In addition, the code takes advantage of compact expressions for 
Higgs boson plus four parton virtual amplitudes that had previously been collected in MCFM. 

The code of our new Higgs boson production generators can be accessed via the POWHEG 
BOX svn repository (see the POWHEG BOX web page http://powhegbox.inib.infn.it for 
instructions). The new MadGraph4 interface is also available there, so that people willing 
to develop their own POWHEG code may benefit from its use. 

A. The interface with MadGraph4: technical details 

The MadGraph4 interface is available under the POWHEG-BOX/MadGraphStuf f directory. To 
use the interface, create a process directory in the POWHEG-BOX directory, and copy the 
whole content of the MadGraphStuf f directory to this new process directory. To specify 
the process, edit the file Cards/proc_card.dat and set the process and the physics model 
as in MadGraph4. Always enter the real emission process, i.e. the Born process plus an 
extra jet j . For example, to generate Higgs boson plus two jets at a proton-proton collider, 
we entered 

P P > H j j j 
QCD=3 
QED=0 
HIG=1 

It is recommended to set the parameters QCD and QED to be exactly equal to the number of 
strong and electroweak interactions in the real-emission process (excluding the couplings 
present in the effective vertex, if any). In the case of the heft model (see below), it is also 
needed to set the parameter HIG=1 to allow for the effective Higgs boson to gluon coupling 
to be included in the Feynman diagrams. 

The ordering of the particles in the process should follow the POWHEG BOX conventions: 

1. first particle: incoming particle with positive rapidity 

2. second particle: incoming particle with negative rapidity 

3. from the third particle onward: final-state particles ordered as follows 

• colourless particles first, 

• massive coloured particles, 

• massless coloured particles. 

The default MadGraph4 interface has been validated to work with the following three physics 
models: 

• sm: this is the default MadGraph4 model with a massive top and bottom quark, 
diagonal CKM matrix and massless electrons and muons. 
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• smckm: this is the same as the sm model, however, the full CKM matrix can be 
specified in the parameter card. 

• heft: this is the same as the default sm model with the inclusion of the effective 
coupling between gluons and the Higgs boson in the large top-quark mass limit. The 
Higgs boson Yukawa interaction with bottom quarks is neglected in this model. 

Although the interface has not been tested with other physics models, it is straightforward 
to extend the interface to work with any simple extension of the Standard Model. 

To generate the process, execute the NewProcess . sh script. This will compile the 
MadGraph4 code, generate the (correlated) Born and real-emission squared amplitudes and 
create the libraries that can be linked to the POWHEG BOX. In particular, libmadgraph. a 
contains the matrix elements, libdhelasS.a the HELAS routines and libmodel.a the 
physics model. Furthermore, the following files are written in the current directory: 

• Born.f : it contains the routines setborn to compute the (correlated) Born squared 
matrix elements and borncolour_lh that assigns the colour flow to a Born process. 

• real . f : it contains the routines setreal to compute the real-emission squared matrix 
elements. 

• Cards/param_card.dat: this is the input card where all the model parameters need 
to be specified. 

• init_couplings . f : this contains the init_couplings routine that sets all the model 
parameters specified in the parain_card.dat. The coupling constants that depend 
on the strong coupling st_alpha or from the event kinematics should be specified in 
the routine set_ebe_couplings, that is updated event-by-event. 

• coupl.inc: it contains the common blocks for all the couplings used by MadGraph4. 

To complete the implementation of a process in the POWHEG BOX, the user must provide 
the Born phase-space in the file Born_phsp . f and the virtual squared matrix elements in 
the file virtual . f . Also, no information on possible intermediate resonances in the matrix 
elements is kept. The user needs to specify explicitly in the routine f inalize_lh (in the 
Born.f file) which resonances should be written in the LHE file, so that the shower Monte 
Carlo program can deal with them correctly. Finally, the resulting code can be compiled 
by executing the command 

$ make pwhg_main 
B. PYTHIA setup 

The sequence of PYTHIA calls we have used in the calculation of the results presented in 
sec. 5 is the following: 

• without hadronization 
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call PYTUNE(320) 

call PYINITCUSER' ," ," ,OdO) 
c Hadronization off 

mstp(lll)=0 
c primordial kt off 

mstp(91)=0 
c No multiple parton interactions 

if (mstp(81) . eq. 1) then 
c Q2 ordered shower 

mstp(81)=0 

elseif (mstp(81) .eq.21) then 
c p_T"2 ordered shower 
mstp(81)=20 

endif 

call PYABEG 
call PYEVNT 
call PYANAL 

• with hadronization 

call PYTUNE(320) 
call PYINITCUSER' , " , " ,OdO) 
c switching off MPI 
mstp(81)=20 
call PYABEG 
call PYEVNT 
call PYANAL 
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